install.packages("sf")Cartography with R - Exercises
Create reproducible maps with R
Vector Data & Geoprocessing
Vector Data
sf package
Install the sf package.
Import and export
Import municipalities.
library(sf)#> Warning in CPL_gdal_init(): GDAL Message 1: GDAL AVIF driver was built against
#> libavif 1.1.1 but is running against 1.2.0. Runtime issues could occur#> Linking to GEOS 3.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE#> WARNING: different compile-time and runtime versions for GEOS found:#> Linked against: 3.13.1-CAPI-1.19.2 compiled against: 3.13.0-CAPI-1.19.0#> It is probably a good idea to reinstall sf (and maybe lwgeom too)st_layers(dsn = "data/lot.gpkg")#> Driver: GPKG
#> Available layers:
#> layer_name geometry_type features fields crs_name
#> 1 communes Multi Polygon 313 12 RGF93 v1 / Lambert-93
#> 2 departements Multi Polygon 96 5 RGF93 v1 / Lambert-93
#> 3 restaurants Point 694 2 RGF93 v1 / Lambert-93
#> 4 elevations Point 5228 1 RGF93 v1 / Lambert-93
#> 5 routes Line String 1054 2 RGF93 v1 / Lambert-93mun <- st_read("data/lot.gpkg", layer = "communes")#> Reading layer `communes' from data source
#> `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 313 features and 12 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
#> Projected CRS: RGF93 v1 / Lambert-93Explore and display
- Plot municpalities with
sf. The filling color must be blue and borders must be thick.
plot(st_geometry(mun), col = "lightblue", lwd = 2)
- Plot municipalities with
mapsf. The filling color must be lightgreen and borders must be white and thin.
# install.packages("mapsf")
library(mapsf)
mf_map(mun, col = "lightgreen", border = "#ffffff", lwd = .5)
- Add the restaurants to the municpalities map.
rest <- st_read("data/lot.gpkg", layer = "restaurants")#> Reading layer `restaurants' from data source
#> `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 694 features and 2 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 542250 ymin: 6350950 xmax: 634198.4 ymax: 6437250
#> Projected CRS: RGF93 v1 / Lambert-93mf_map(mun, col = "lightgreen", border = "#ffffff", lwd = .5)
mf_map(rest, pch = 21, col = "red", cex = 1, add = TRUE)
- How many municipalities are there?
# dim(mun)[1]
nrow(mun)#> [1] 313# length(unique(mun$INSEE_COM))- Which commune has the largest population?
# mun[order(mun$POPULATION, decreasing = T), 2][1, ]
mun[mun$POPULATION == max(mun$POPULATION), 2]#> Simple feature collection with 1 feature and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 570470.4 ymin: 6368056 xmax: 581068.4 ymax: 6380468
#> Projected CRS: RGF93 v1 / Lambert-93
#> NOM_COM geom
#> 39 Cahors MULTIPOLYGON (((580994.5 63...Attribute selection and join
- Import the com.csv file.
mun_df <- read.csv("data/com.csv")- Join this dataset to the municipality layer.
mun_final <- merge(
x = mun, # the sf object
y = mun_df, # the data.frame
by = "INSEE_COM", # identifier in x and y
all.x = TRUE # keep all lines
)- Extract municipalities with more than 500 active workers and whose share of active workers in the total population is greater than 30%.
mun_sel <- mun_final[mun_final$ACT > 500 & mun_final$SACT > 30, ]- Plot a map with all municipalities in gray and selected municipalities in red.
mf_map(mun_final, col = "grey40", border = 'lightblue')
mf_map(mun_sel, col = "red", add = TRUE)
Geoprocessing
Spatial selection and join
- Import the layers of municipalities and restaurants.
rest <- st_read("data/lot.gpkg", layer = "restaurants")#> Reading layer `restaurants' from data source
#> `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 694 features and 2 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 542250 ymin: 6350950 xmax: 634198.4 ymax: 6437250
#> Projected CRS: RGF93 v1 / Lambert-93mun <- st_read("data/lot.gpkg", layer = "communes")#> Reading layer `communes' from data source
#> `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 313 features and 12 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
#> Projected CRS: RGF93 v1 / Lambert-93- Perform a spatial join to find the name and identifier of the municipality where each restaurant is located.
resto <- st_join(x = rest, y = mun[, "INSEE_COM"],
join = st_intersects,
left = FALSE)
resto#> Simple feature collection with 694 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 542250 ymin: 6350950 xmax: 634198.4 ymax: 6437250
#> Projected CRS: RGF93 v1 / Lambert-93
#> First 10 features:
#> id nom INSEE_COM geom
#> 1 212 forest_cobra 46154 POINT (611092.9 6365662)
#> 2 2327 meteoropathologic_marlin 46003 POINT (596050 6414450)
#> 3 2328 pumice_unicorn 46003 POINT (596350 6414950)
#> 4 2329 indifferent_nakedmolerat 46003 POINT (595750 6417150)
#> 5 2330 pedagogic_eskimodog 46004 POINT (613350 6405150)
#> 6 2331 painstaking_alaskanhusky 46005 POINT (556763.5 6377394)
#> 7 2332 taxonomical_archaeopteryx 46006 POINT (573350 6411050)
#> 8 2333 spotless_hornedtoad 46008 POINT (560650 6389150)
#> 9 2334 dramatic_hyrax 46009 POINT (610450 6398050)
#> 10 2335 carefree_mountainlion 46009 POINT (610640.9 6397839)Geometrical operations
- Compute the number of restaurants per municipality.
mun$n <- lengths(st_intersects(mun, resto, sparse = T))
########### OR
# resto <- st_join(x = rest, y = mun[, "INSEE_COM"],
# join = st_intersects,
# left = FALSE)
# rest_by_mun <- aggregate(x = list(n = resto$id),
# by = list(INSEE_COM = resto$INSEE_COM),
# FUN = length)
# mun <- merge(mun, rest_by_mun, "INSEE_COM", all.x = TRUE)
# mun[is.na(mun$n), "n"] <- 0- Which municipalities have more than 10 restaurants and fewer than 1000 inhabitants?
mun_sel <- mun[mun$n > 10 & mun$POPULATION < 1000, ]- Create a map showing all the municipalities in grey and the municipalities selected above in red.
mf_map(mun_final, col = "grey40", border = 'lightblue')
mf_map(mun_sel, col = "red", add = TRUE)
Cartography
The Cartographic Language
How to represent the following variables:
- A municipal population
Proportional Symbols - A median age by region
Choropleth map, ordered colors
- A population growth rate
Choropleth map, ordered colors
- The administrative status of municipalities
Typologie map, unordered colors
- Life expectancy per country
Choropleth map, ordered colors
mapsf
Map Types
- Import the com.csv file.
mun_df <- read.csv("data/com.csv")- Join this dataset to the municipality layer.
mun <- st_read("data/lot.gpkg", layer = "communes")#> Reading layer `communes' from data source
#> `/home/tim/Documents/prj/cartography_with_r/data/lot.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 313 features and 12 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
#> Projected CRS: RGF93 v1 / Lambert-93mun <- merge(
x = mun, # the sf object
y = mun_df, # the data.frame
by = "INSEE_COM", # identifier in x and y
all.x = TRUE # keep all lines
)- Create a map of the working population.
library(mapsf)
mf_map(mun, col = "ivory4", border = "ivory" , lwd = .3)
mf_map(
x = mun ,
var = "ACT",
type = "prop",
inches = .35,
symbol = "square",
border = "white",
lwd = .75,
col = "#881094",
leg_title = "Working Population",
leg_pos = c(622700, 6372800)
)
mf_title("Working Population in the Lot Region")
- Create a map of the share of the working population in the total population.
mf_distr(mun$SACT)
summary(mun$SACT)#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 7.812 23.104 26.754 27.265 31.537 52.083# Classification method based mean and standard deviation
bks <- mf_get_breaks(mun$SACT, breaks = "msd", central = FALSE)
# Double color palette
cols <- mf_get_pal(c(3,4), palette = c("Mint", "Burg"))
mf_map(
x = mun,
var = "SACT",
type = "choro",
breaks = bks,
pal = cols,
border = "ivory4",
lwd = .2,
leg_title = "%",
leg_val_rnd = 0
)
mf_title("Share of the working population in the total poulation")
Map Layout
- Create a map showing the working population in the industrial sector.
- Add the necessary layout elements.
- Make the map more intelligible, more explicit.
- Export the map in PNG format, 1000 pixels wide.
Workers + Workers in the industrial sector
dep <- st_read("data/lot.gpkg", layer = "departements", quiet = TRUE)
# Workers + Workers in the industrial sector
par(mfrow = c(1, 2))
mf_map(x = mun, border = "white", lwd = .2, add = F)
mf_map(x = dep, lwd = 1, col = NA, add = TRUE, lend = 0)
mf_map(x = mun, var = "ACT",
type = "prop",
inches=.2,
leg_frame = T,
leg_title = "N.")
mf_title("Working Population", tab = FALSE)
mf_credits(
paste0("Admin Express COG Carto 3.0, IGN - 2021 & ",
"BD CARTO® 4.0, IGN - 2021 ; Recensements harmonisés - ",
"Séries départementales et munmunales, INSEE - 2020\n",
"Auteurs : T. Giraud, 2025"), bg ="#ffffff80")
mf_map(x = mun, border = "white", lwd = .2, add = F)
mf_map(x = dep, lwd = 1, col = NA, add = TRUE, lend = 0)
mf_map(x = mun, var = "IND",
type = "prop",
inches=.2,
leg_frame = T,
leg_title = "N.")#> 50 '0' values are not plotted on the map.mf_title("Working Population in the Industrial Sector", tab = FALSE)
mf_scale(5)
mf_arrow(pos = "topright")
Share of workers in the industrial secto
mf_distr(mun$SACT_IND)
summary(mun$SACT_IND)#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.000 7.407 14.286 16.282 23.077 100.0002 municipalities have a 100% share of workers in the industrial sector. Both have less than 15 workers. They are outliers.
# Remove outliers
mun_sel <- mun[mun$ACT > 15, ]
mf_distr(mun_sel$SACT_IND)
summary(mun_sel$SACT_IND)#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.000 7.408 14.286 15.783 22.466 60.000bks <- mf_get_breaks(mun_sel$SACT_IND, nbreaks = 5, breaks = "quantile")
mf_map(x = mun,
var = "SACT_IND",
type = "choro",
breaks = bks,
leg_val_rnd = 0,
pal = "Red-Yellow",
leg_title = "%",
add = FALSE,
col_na = "grey",
leg_no_data = "less than 15 workers")#> 2 values are outside the class limits and will be classified as 'NA'.mf_title("Share of workers in the industrial sector")
mf_scale(5)
mf_arrow(pos = "topright")
mf_credits(paste0("Auteurs : T. Giraud, 2025\n","Admin Express COG Carto 3.0, IGN - 2021 & ",
"BD CARTO® 4.0, IGN - 2021 ;\nRecensements harmonisés - ",
"Séries départementales et munmunales, INSEE - 2020"))
Final Map
mf_export(x = mun, filename = "img/map.png", width = 1000, res= 120)
mf_map(x = mun, border = "white", lwd = .2, add = T)
mf_map(x = dep, lwd = 1, col = NA, add = TRUE, lend = 0)
mf_map(mun, c("ACT", "SACT_IND"), "prop_choro",
breaks = bks,
pal = "Red-Yellow",
inches = .4,
border = "white", lwd = .7,
leg_val_rnd = c(0,1),
leg_pos = c(538000,6442000, 538000, 6424000),
leg_title = c("N. workers",
"Workers in the\nIndustrial Sector\n(%)"),
col_na = "grey",
leg_no_data = "Less than 15 workers")#> 2 values are outside the class limits and will be classified as 'NA'.# Ajout d'annotations
mf_annotation(x = mun[mun$NOM_COM=="Biars-sur-Cère",], txt = "Andros\n(jam factory)",
halo = TRUE, cex = 1)
mf_annotation(x = mun[mun$NOM_COM=="Figeac",], txt = "Aerospace\nIndustry",
pos = "bottomright", halo = TRUE, cex = 1)
mf_annotation(x = mun[mun$NOM_COM=="Gramat",], txt = "La Quercynoise\n(duck meat factory)",
pos = "topleft", halo = TRUE, cex = 1)
mf_title("MECANIC VALLEE")
# inset->
mf_inset_on(fig = c(.8,0.98,0.1,0.3))
mf_map(dep, lwd = .1)
mf_map(mun, border = NA, add = TRUE, col = "#D7FF68")
box(lwd = .5)
mf_inset_off()
# <- inset
mf_scale(5)
mf_arrow("topright")
mf_credits(paste0("Auteurs : T. Giraud, 2025\n","Admin Express COG Carto 3.0, IGN - 2021 & ",
"BD CARTO® 4.0, IGN - 2021 ;\nRecensements harmonisés - ",
"Séries départementales et munmunales, INSEE - 2020"))
dev.off()#> png
#> 2